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Non-Gaussianity in the distribution of inflationary perturbations, measurable in statistics of the 
cosmic microwave background (CMB) and large scale structure fluctuations, can be used to probe 
non-trivial initial quantum states for these perturbations. The bispectrum shapes predicted for 
generic non-Bunch-Davies initial states are non-factorizable (“non-separable”) and are highly oscil¬ 
latory functions of the three constituent wavenumbers. This can make the computation of CMB 
bispectra, in particular, computationally intractable. To efficiently compare with CMB data one 
needs to construct a separable template that has a signihcant similarity with the actual shape 
in momentum space. In this paper we consider a variety of inflationary scenarios, with different 
non-standard initial conditions, and how best to construct viable template matches. In addition 
to implementing commonly used separable polynomial and Fourier bases, we introduce a basis of 
localized piecewise spline functions. The spline basis is naturally nearly orthogonal, making it easy 
to implement and to extend to many modes. We show that, in comparison to existing techniques, 
the spline basis can provide better hts to the true bispectrum, as measured by the cosine between 
shapes, for sectors of the theory space of general initial states. As such, it offers a useful approach 
to investigate non-trivial features generated by fundamental properties of the inflationary Universe. 


I. INTRODUCTION 


We are fortunate to live in a time when cosmologi¬ 
cal datasets can probe the Universe in exquisite detail. 
In particular, the cosmic microwave background (CMB) 
provides a rich source of information about the very early 
Universe, and is an especially precise probe of the infla¬ 
tionary paradigm. An important question that we are 
now in a position to probe, more thoroughly than ever, 
is: what is the initial quantum state of inflationary fluc¬ 
tuations? It is usually taken to be the Bunch-Davies 
state, but from the point of view of treating inflation 
as an effective theory it is not unreasonable to consider 
the choice of state to be open, subject to the conditions 
that it allow for inflation to occur and that it be con¬ 
sistent with field theoretic precepts. Explicit examples 
of scenarios that give rise to non-Bunch-Davies initial 
conditions for inflation can be found in m\- Assum¬ 
ing that the initial state is more general than the free 
vacuum, such as a Bogoliubov transform of the Bunch- 
Davies state or even a mixed state, we can calculate its 
imprint on cosmological observables like the power spec¬ 
trum and bispectrum of inflationary perturbations, and 
in turn those of CMB temperature anisotropies EmEi- 
122] • Whether these effects can actually be observed in 
cosmological data depends on the extent of departure 
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from a Bunch-Davies state, the number of e-folds of in¬ 
flation beyond the minimum required, and of course the 
sensitivity of our experiments |23| . 

Given these choices, how can we narrow down the pos¬ 
sibilities? The initial state of perturbations will leave 
its imprints on their correlation functions, for example, 
logarithmic oscillatory modulations in the CMB power 
spectrum [21]. Higher order correlators, such as the bis¬ 
pectrum p5l |2S|, are extremely sensitive to deviations 
from the Bunch-Davies state. The bispectrum carries in¬ 
formation both about the amplitude of the correlations, 
typically encoded in /nlj as well as of preferred config¬ 
urations in momentum space; the three momenta must 
form a triangle, but the shape of the triangle is sensi¬ 
tive to both the interactions of the inflaton as well as the 
initial state. 

The bispectrum for general initial states is highly os¬ 
cillatory and cannot be written in a separable form, i.e. 
as a product of separate functions of the three momen¬ 
tum modes. This makes the study of these states via the 
bispectrum computationally difficult. For such shapes 
we usually construct a basis of separable functions and 
rewrite the desired shape as a sum over many such ba¬ 
sis functions, for example, using a polynomial basis m, 
Fourier basis [2H], or divergent basis As long as 

the original and reconstructed shapes are very similar 


^ In special shape-specific cases, other basis sets can be used: 
for example, feature or resonant models exhibiting linear or 
logarithmic oscillations, respectively, can be efficiently recon¬ 
structed through a one-dimensional (ID) expansion in the sum 
of wavenumbers, fci -j- /c 2 + [301EJ. 
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and have a significant overlap, or in other words have a 
cosine close to unity, we can look for signatures of the 
reconstructed shape in the CMB bispectrum and be as¬ 
sured that the result will be an accurate reflection of what 
we would have obtained for the real shape. Non-Bunch- 
Davies shapes, however, can be difficult to efficiently de¬ 
scribe with existing bases |26fE5] . This leaves open the 
possibility that signatures of general initial states could 
still be present, yet undetected, in the CMB data. 

In this paper we describe a new basis of piecewise spline 
functions and use it to fit non-Bunch-Davies shapes of 
the bispectrum. The spline basis consists of polynomial 
functions defined locally, between various points called 
“knots” in 3D space |32l ESI • This makes it particularly 
suitable for smoothing and interpolating data with com¬ 
plex patterns, and in our case, for capturing localized 
features of any complicated shape of the bispectrum. An¬ 
other immediate advantage of defining localized functions 
is that the basis functions are orthogonal to a good ap¬ 
proximation, and there is no need to perform a Gram- 
Schmidt-like orthogonalization process on the basis func¬ 
tions. This makes the spline basis easy to implement and 
to extend, using a large number of mode functions to 
capture fine features in the bispectrum. We find that the 
spline basis performs at least as well as the polynomial 
and Fourier bases in describing most non-Bunch-Davies 
shapes, and for many shapes offers significant improve¬ 
ments (in the cosines) over existing techniques. 

The remaining paper is organized as follows. In section 
[IT] we briefly review how we define the initial conditions 
for the perturbations and obtain the bispectrum for gen¬ 
eral initial states. We describe the spline basis in section 


cussion on the scope of this work. Appendices and 
11 contain details on calculations of the correlation func¬ 
tions for general initial states, and appendix [C] describes 
our numerical implementation of the spline basis. 


III and use it to fit non-Bunch-Davies shapes in section 

IV We conclude in section with a summary and dis- 


II. NON-BUNCH-DAVIES SHAPES 

We usually describe primordial correlation functions 
in terms of the curvature perturbation since this 

quantity is conserved outside the horizon [33|. It is de¬ 
fined as the perturbation in the local scale factor a{t); 
the metric perturbation hij {t, x) is then written as hij = 
a?e^^5ij. In an effective field theory setting, C(C^) is re¬ 
lated to the Goldstone mode of time reparameterization 
symmetry breaking, usually denoted as 7r(t,ai) |351 136|. 
In this section we describe how the choice of initial state 
for the perturbation (^{t,x) (or equivalently Tr(t,x)) af¬ 
fects the primordial bispectrum fc 2 , defined via 

(CfciCfc2Cfc3) = {2n)^6^{ki + k2 + k3)Bt:{ki,k2,h). 

Starting with the Einstein-Hilbert action, we can cal¬ 
culate the action for scalar perturbations directly in 
the C“g 8 -uge. Writing down the most general Lorentz- 
invariant scalar-tensor theory with second order equa¬ 


tions of motion results in the Horndeski action for the 
perturbations For P{X, (p) models of inflation, 

with X = —g^^d^(l)d''(l>,^ the leading order in slow-roll 
Horndeski action at cubic order in the perturbations is 
given by 


S = 


d^xdt a^< 




VAiC^ + AaCC" 



^boundary • 


( 1 ) 


Here we have set Mpi = 1, e = —H/H^ is the slow-roll 
parameter, Cg is the effective sound speed for perturba¬ 
tions, and the couplings Ai — A 3 are given by 


Ai = 

Hct[ * E j ’ 

(2) 

A 2 = 


(3) 

A 3 = 


(4) 


with A = andE = = 

t_x2 

. The boundary terms in eq. (1) are important and in 
general do contribute to the bispectrum when the initial 
state is different from the Bunch-Davies vacuum; here we 
assume for simplicity that the initial state for C{t,x) is de¬ 
fined in such a way as to cancel these boundary terms. It 
is also worth noting that the action in eq. is equivalent 
to that obtained using an effective field theory approach 
[351136| up to boundary terms, and the bispectrum for 
general initial states agrees between the two actions, as 
shown in |42| . 

The usual method to obtain correlation functions of 
C(t, x) is to use in-in perturbation theory |43f[5D] . The 
initial state can be input as a density matrix defined at 
the initial time to | 20 | ; for simplicity we will assume that 
the initial density matrix is pure and Gaussian.^ Dou¬ 
bling the fields on the plus and minus branches of the 
in-in contour, we can calculate the Green’s function cor¬ 
responding to the quadratic (or “free”) part of the action 
in eq. 0 - including the effect of a general initial state; 


^ The Horndeski action does not describe ghost inflation, however, 
which can be included in an effective field theory setting. An¬ 
other example outside the Horndeski domain is Horava-Lifshitz 
gravity, in which Lorentz invariance is explicitly broken. 

^ The action for many single scalar field models of inflation can be 
written as 5" = d'^Xy/—g [/? + 2P{X, 0)] , with (p controlling 
the dynamics of both the background and perturbations. 

^ Here by “pure” we are distinguishing between pure and mixed 
quantum states. Mathematically speaking, Tr(p^) = 1 for pure 
states, while Tr (p^) < 1 for mixed states, p being the density 
matrix. Relaxing the pure state assumption leads to qualitatively 
very similar results to what we discuss here. By “Gaussian” we 
mean that the action describing the initial state is quadratic. 
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details of this calculation can be found in the appendix 
of (2D]. The Green’s function for C{t,x) is found to be 




Cs 1 { G++{r],7j') G+ 

2e a{T])a{r]') \ J ’ 

(5) 


where rj is the conformal time defined as ij = f dt/a and 
the factor out front comes from rewriting the quadratic 
action in terms of the canonically rescaled field 
C = :^^X- The functions G^’^{ri,rj') are given by 


Gt^{rj,rj') 

= fkiri)fkivy{ri-v') 



+ fkiv)fkiv'W{v'-v), 

( 6 ) 

G+-{rj,rj') 

= fkiv)fkiv'), 

(7) 

Gk^iv,v') 

= G^*(77 , 77'), 

( 8 ) 

Gk~iv.v') 

= G++*(77, 770, 

(9) 


and the mode functions f^iv) are solutions to the sec¬ 
ond order differential equation resulting from the Green’s 
function equation, 


fkiv) = a^hkirj) + P^hlirt) 

- ‘2fkivo)Akgk{v)0{vo - v) ■ ( 10 ) 


The last term above is an additional contribution from 
the initial density matrix, Ak being the kernel that mul¬ 
tiplies the term in the initial state action HOI; it 

does not, however, contribute to the bispectrum, and so 
we will ignore it in our discussion. The Bogoliubov co¬ 
efficients a^, are functions of kernels in the initial 
density matrix, the mode functions, and time derivatives 
of the mode functions, all at the initial time, and satisfy 
«fc = ^k = 1^ - l/3fe 1^ = 1- Finally, the 

function hk{rj) is defined at leading order in slow-roll as 

hkiv) = , (11) 


both choices of initial time in section |IV| when we ap¬ 
ply the spline basis to non-Bunch-Davies shapes of the 
bispectrum. 

Let us now discuss how general initial states modify the 
bispectrum. Observables, such as the bispectrum, can be 
calculated using any combination of plus and minus fields 
on the in-in contour. For Gaussian initial states, we can 
calculate the three-point function in the perturbations as 


i 


X exp 


,( 12 ) 


where is the cubic part of the action in eq. Q writ¬ 
ten in momentum space, with conformal time derivatives, 
and with the time integral running from rjo to rj. The 
subscript “G” indicates that Wick contractions on the 
right are carried out using the Gaussian theory. At lead¬ 
ing order, only the three operators C^, and C(9iC)^ 
contribute, and we can write the three-point function at 


late times ^ ^ contribu¬ 

tions from each of these three operators. To calculate the 
three-point function in eq. at late times we need the 
function G^.’’'’^( 0 , 77 ') and its derivative 9,,'G^’^''’(0, 77 ') 
for 77 ' > 770 ; using eqs. (discarding the 0(770 — 77 ) 

term), and ( 11 ) these are given by 


Gi’++(0,77') = 




-I- 6fc(l -I- 7Csfc77')e“*°='"’' 


(13) 


and 


9^,Gi’++(0,77') = 


^ ^^icskrj' _|_ ^^^—icskri' 


Aek 


Qj}iC 




(14) 


where we have defined the functions 


where .^ 3/2 i® ^ Hankel function. The Bunch-Davies 
choice consists of setting the initial time 770 —>■ — 00 , the 
initial density matrix to unity, and additionally = 1 , 
= 0 so that the mode function fjP (rj) picks out the 
positive frequency solution proportional to 

The time 770 at which the initial conditions are set can 
be taken to be a constant time in the past at the on¬ 
set of infation, or can be considered as a scale-dependent 
quantity, 770 (fc). In the latter case, the initial conditions 
for each k mode are set at the time when the physical 
momentum corresponding to this mode Csk/a{r]o) (with 
^ivo) = ~^/iVoH) during inflation, at leading order) 
crosses a fixed energy scale A of new physics. The Bo¬ 
goliubov transform in eq. (101 is correct for either choice 
of TjQ. In appendix]^ we show that for a scale-dependent 
initial time, this solution leads to the well-known oscil¬ 
lations in the late-time power spectrum uniin]. We use 


ak = (ttfe - Pk ) , 

h = - 


(15) 

(16) 


Using these in eq. (12 1 and performing the time integrals 


we find that the contributions to the bispectrum from 
the three operators are given by 




32 


(2^)3 <53 


X Ai 




/ . ‘'fci % 

/,m,n =0 


X J-^3 ((-l)'fci,(-l)"^fc2,(-l)”fc3,77o) 


-I- c.c., 


(17) 
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32 


(2Trf fc,) 


X A 2 




1 



_ 1,771,71—0 



C 


(n) 

^3 


X ((-iyk,,(-irk2,(-irks,r/o) 

+ c.c. , (18) 


and 


(n) 

\^fei^fe2^fc3/c(ac)" 


I 7?—V 0 “ 


1 




X A. 




E 

Z,?7i,n=0 


(0 J' 


ml {n 1 


Cl cl 


X ((_!)'*,, (_i)™fc2,(-l)”A3,%) 

+ c.c., (19) 

where “c.c.” denotes complex conjugate and 


Ad 


Ofe 1 = 0 

bk i = 1. 


( 20 ) 


The functions Jy/ 2 , and are written out 

explicitly in appendix]^ The above equations give us 
the leading order result for the bispectrum for general 
initial states. The non-Bunch-Davies contributions to 
the bispectrum are strongly peaked in the flattened limit 
/ci « ^2 + ks (assuming that ki is the largest momentum 
mode) and in the squeezed limit k^ ^ ki k, k^- (We show 
these enhancements and discuss apparent divergences in 
both limits in appendix i) Further, these shapes are 
highly oscillatory, which makes them even harder to con¬ 
strain using the CMB bispectrum. In the next section 
we discuss the spline basis that we use to rewrite these 
shapes as a sum of separable functions. 


III. THE SPLINE BASIS 

B-splines, short for “basis splines”, are a well- 
established, and conceptually simple, mathematical for¬ 
malism for curve fitting, using a set of piecewise poly¬ 
nomial functions [321 [33]. For a basis in one dimension, 
one chooses a set of “knots”, {xq, xi,..., xn}, represent¬ 
ing the points at which the polynomial function pieces 
are joined, and the degree q of the polynomials. For ex¬ 
ample, a ID spline basis spanning 0 < x < 1 with a set 
of sbc piecewise cubic polynomials is shown in fig. As 
an explicit example of the functional form of the basis, 
one of the basis functions shown in the figure is 

r I (4x - 18x2 -F 21x3) 0 < X < i 

^ I i (8 — 36x -I- 54x2 — 27x3) | < x < |. 

( 21 ) 



FIG. 1: Example spline basis generated from knots at x = 
{0, 0, 0,0,1/3, 2/3,1,1,1,1} and polynomials of degree q — 3. 


The numerical coefficients defining the spline basis for in¬ 
put knots and q are easily generated using existing soft¬ 
ware libraries in many languages. 

A general ID function, /(x), can be expanded, and 
approximated, using the spline basis. 


N-l 

/'(x) = ^ a„B„(x). (22) 

n—0 

The expansion coefficients are computed by first tabu¬ 
lating a sample of (xi, fi) values, where fi = f{xi), cor¬ 
responding to M data points, and finding the values of 
{an} that minimize the least-squares function, 

M / N-i \ 2 

LS = E - E ^nBnixy . (23) 

2—1 \ n —0 / 

In practice, this requires solving the linear system of 
equations given by Ba = B"^ f, where B here is an 
M X N matrix containing the values of Bn{xi). 

Similarly, we can generate a higher-dimensional basis 
by multiplying ID b-splines together. For example, a 3D 
b-spline fit for a shape® can be expressed as 

— 1 k j 

s'ik.MM = eee Aijk[Byki)Bj{k2)Bkik3) 

k—0 j—0 2—0 

-I- perms] , (24) 

where “perms” refers to permutations of fci, fc 2 , and k^. 
We note that the number of relevant {i,j,k) combina¬ 
tions in such a fit is not N^, because each spline mode is 
symmetric in the three wavenumbers, and not all (z, j, k) 
combinations will correspond to a spline mode that is 


® A shape function is usually denoted by 5(A:i, fc 2 , fcs), not to be 
confused with the action S. 
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non-zero for wavenumber combinations that can form a 
closed triangle. In the 2D and 3D cases, efficient algo¬ 
rithms solving for the expansion coefficients already 
exist in the literature m- We illustrate how to use the 
spline basis further through explicit examples and include 
snippets of our numerical codes in appendix [C| 

In contrast to other bases used to create non-Gaussian 
templates made of globally varying functions, such as 
polynomial, Fourier, and divergent functions, the spline 
basis consists of a set of localized modes, each of which 
describes only a small region of the allowed fc-space. 
This makes the spline basis very well-suited to describ¬ 
ing shapes such as those of non-Bunch-Davies models, 
that are characterized by highly-peaked features concen¬ 
trated on very flattened triangles. Because of the local¬ 
ization of spline modes, the modes are by construction 
nearly orthogonal, and no orthogonalization procedure 
(such as Gram-Schmidt) is used in our implementation. 
Avoiding any explicit orthogonalization is an advantage, 
as orthogonalizing the polynomial/Fourier basis sets via 
Gram-Schmidt is numerically unstable, and requires very 
high precision throughout to create a basis with a large 
number of modes.® 


IV. FITTING NON-BUNCH-DAVIES SHAPES 

In this section we consider two non-Bunch-Davies 
shapes included in Planck’s analysis that were not well- 
reconstructed in multipole-space using the polynomial 
modal expansion, S'nbdi and S'nbd2 We also con¬ 

sider a set of generalized non-Bunch-Davies shapes, S^s, 
S(-^ 2 , and with different assumptions for the ini¬ 

tial time boundary CsTJo and its potential wavenumber de¬ 
pendence, Cs|77o| = 10^ Mpc, {A/H)/{ki + k 2 + fcs), and 
{A/H)/ki, where we allow A/H to be 10 or 10®, rep¬ 
resenting a physically-motivated range of values. These 
shapes span a wide variety of possible non-Bunch-Davies 
features in the bispectrum, and allow us to compare our 
results with previous analyses. 

For any non-separable primordial shape S, and given a 
choice of basis {M„}, one can compute a separable fit S' 


® The numerical instability of the classical Gram-Schmidt algo¬ 
rithm can be partially mitigated by instead adopting the mod¬ 
ified Gram-Schmidt algorithm, which is what we have imple¬ 
mented in generating the polynomial and Fourier modes. 

^ In the notation of 1261 . Sned* here is equal to 

(fclfc2A:3)2BNBDY i = 1,2. 

® The shapes and S'^(a ^)2 here are related to 

corresponding leading order non-Bunch-Davies correc¬ 
tions to the bispectrum; for example, (fci, ^ 2 , ^ 3 ) = 

(fcl/C2^3)^ + -^<^3 (^1,-^2, ^3,’?0) + 

(All,/c2) —^3) ??o) + c.c.j, where we have set = 1 

and bfc- = 0.01. In general b^. can have some scale-dependence 
as long as it does not spoil constraints from backreaction of the 
energy density in the initial state. 


that approximates (S' as a iinear combination of separable 
basis functions, 


S"(fci,fc2,fc3) = ^a„7W„(/ci,A:2,fc3). (25) 

n 


The similarity between the original shape and its fit is 
quantified by a cosine. 


cos{S, S') 


{S,S') 

y'{S,S) {S', S') ’ 


(26) 


where the inner product is defined in Fourier space with 
a choice of weighting, 

{S,S') = f dVr S{k,,k2,k3)S'{kuk2,k3)w{k,,k2,k3) . 

(27) 


The volume Vt includes only those combinations of ki, 
k 2 , and k^ that can form a closed triangle, with each 
wavenumber satisfying fcmin < fci, ^ 2 , fca < ^maxi where 
fcmin = 10“® Mpc“® and fcmax = 0.1 Mpc“^. The 
weight w{ki,k 2 ,k^) is typically taken to be either unity 
or l/(A:i + k 2 + k^), where the latter choice is meant 
to represent a more accurate reflection of the scaling of 
the covariance of the GMB bispectrum, such that the 
Fourier-space cosine is closer to the multipole-space co¬ 
sine between the GMB bispectra corresponding to the 
shapes S and S'. However, in general we find that both 
choices result in similar cosines. In our analysis, we use 
a unit weight for the 3D Fourier basis fits and all 2D fits, 
and a weight of l/(fci - 1-^2 + k^) for the remaining fits. 

In this section we implement the existing polynomial 
and Fourier basis methods, and additionally our new ba¬ 
sis of piecewise splines, to obtain separable fits to non- 
Bunch-Davies shapes. In each case, we quantify the per¬ 
formance of the basis by computing the cosine as a func¬ 
tion of an increasing number of modes used in the fit. 

We use the polynomial basis described in m and the 
Fourier basis described in [^. In each case, three ID 
functions based on either polynomials or sines/cosines of 
ki are multiplied together to form 3D separable functions, 
that are then orthogonalized using a Gram-Schmidt algo¬ 
rithm to create a basis of 3D orthonormalized and separa¬ 
ble functions, called {77,„} for polynomials and {J>i} for 
Fourier modes. Since the basis functions are orthonormal 
in each case, the expansion coefficients {«„} can be com¬ 
puted through inner products between S and the basis 
functions, a„ = (S', 77,„) or a„ = {S,Tn)- 

Alternatively, in the spline basis expansion, three ID 
piecewise spline functions are multiplied together to form 
a basis of 3D separable spline functions {Bn}. In this 
case, each mode is highly localized in a region of Fourier 
space, so any two modes are orthogonal by construction 
unless they have peaks that overlap. While the lack of 
strict orthogonality means that the expansion coefficients 
in the spline basis cannot be computed using simple inner 
products, existing algorithms can solve for the coefficients 
efficiently m- 
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In the subsections that follow, we present polynomial, 
Fourier, and spline basis fits to a variety of non-Bunch- 
Davies shapes. 

A. 5'nbdi and S'nbd2 

The cosines for polynomial, Fourier, and spline fits to 
these shapes are shown in the upper left panel of fig. 
1^ We find that the polynomial expansion produces a 
better fit than the Fourier basis, and polynomial cosines 
typically increase slowly beyond about 100 modes for 
S'nbdi and 50 modes for 5'nbd2) indicating that lower 
order modes contribute most to the fits. While it is pos¬ 
sible that increasing the number of polynomial or Fourier 
modes will increase the cosines further, during our analy¬ 
sis we found that generating large polynomial and Fourier 
basis sets is computationally very demanding. The sep¬ 
arable modes are orthogonalized using a Gram-Schmidt 
algorithm, which is known to be numerically unstable, 
and the instabilities become more severe as higher order 
polynomials are used. On the other hand, the spline ba¬ 
sis expansion as we have implemented it does not require 
orthogonalization, so in comparison a large number of 
modes can easily be generated and used in the separable 
fits. In the case of 5'nbdi and 5nbd2) we find that the 
spline expansion performs similarly well as the polyno¬ 
mial expansion when 200 modes are used, and better fits 
can be achieved by using a larger number of modes. 

B. 5^3, and with fc-independent Csr]o 

These shapes have oscillatory features that grow with 
the size of the triangle (i.e. with ki) and in the flattened 
and squeezed limits. With Cs|77o| = 10^ Mpc, the set of 
200 polynomial modes produces a cosine of 0.52 for 5^3 
and 0.80 for 5^^2, as shown in the upper right panel of 
fig-i" While these modes produce a higher cosine than 
a set of 200 spline functions, better fits can be generated 
with a larger set of splines. 

C. S^ 3 , 5^^2, and 5(‘(a^)2 with fc-dependent CsTIo 

For small A/H, we find that the shapes are easily re¬ 
constructed with both the polynomial and spline meth¬ 
ods using < 200 modes, as shown in lower left panel 
of fig. [I The oscillatory features are of low enough fre¬ 
quency that lower order polynomials are sufficient to cap¬ 
ture most of the features of these shapes, and the spline 
functions also do not need to have a very fine resolu¬ 
tion. If the initial conditions are set at 10//ci, rather 


® The shape for is sufficiently similar to that its fits 

are not shown separately in the figures. 


than 10/(fci + k 2 + k^), then the oscillatory features are 
of somewhat higher frequency; for either shape, however, 
the large global features still allow both the polynomial 
and spline bases to efficiently reconstruct these shapes. 

For cases with fc-dependent CsTJo and larger values 
of A/H, the oscillations have a much higher frequency 
than what can be captured by polynomials, and we 
find that the spline reconstructions perform similarly 
poorly. However, the /c-dependence of CsTJo, whether it 
be (A/i?)/(fci-I-/c2-f fcs) or {A/H)/ki, allows the shapes 
to be rewritten as functions of only two free parameters: 
the ratios x = k^/ki and y = k^/ki. In terms of x and 
y, the oscillation frequency does not increase drastically 
throughout the allowed parameter space, which makes it 
easier to generate good 2D fits. The advantage of rewrit¬ 
ing scale-invariant shapes such as these in terms of two 
free parameters for computing fc-space cosines and CMB 
bispectra has been discussed in earlier works [52] . 

We generate the 2D fits by defining a 2D analogue of 
the polynomial/spline basis sets and the cosine. We find, 
however, that to generate the same number of polynomial 
modes as in the 3D case (200) requires using higher order 
ID polynomials, which exacerbates the numerical issues 
that we encountered in the 3D case, so for the 2D fits we 
only use 100 polynomial modes with unit weight. We do 
not encounter similar numerical issues in the spline fits. 

We show the results for the 2D fits in the Cs|77o| = 
10^/(fci -I- ^2 -I- fca) case in the lower right panel of fig. 
for Cs|??o| = 10^/fci the shapes have very similar features, 
with the oscillations being slightly more rapid and more 
difficult to represent in the latter case. 2D spline fits 
achieve similar cosines as 2D polynomial fits using the 
same number of modes 100), while also being able 
to produce higher cosines through the addition of more 
modes without running into numerical issues, as shown 
in fig. 

V. DISCUSSION 

If the CMB is to be used to its full potential to eluci¬ 
date the details of inflation, we need to be able to both 
effectively characterize a wide variety of potential infla¬ 
tionary signatures while also ensuring that this informa¬ 
tion is accessed in a timely manner. For non-Gaussian 
signatures, present in the CMB bispectrum, a potential 
bottleneck in this process is the production of separable 
templates that accurately match the characteristic prop¬ 
erties of the underlying shape. This issue is particularly 
acute for initial states that differ from the Bunch-Davies 
state, which can generically have highly oscillatory fea¬ 
tures whose characteristic scale can vary over the bispec¬ 
trum configuration space. 

In this work we have analyzed a variety of functions to 
generate templates for bispectrum shapes arising in a va¬ 
riety of non-Bunch-Davies scenarios. We have quantified 
how well different choices of separable basis functions, de¬ 
rived from polynomials, Fourier functions, and b-splines. 
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FIG. 2: Comparisons of shape reconstructions for 5'nbdi (solid) and Snbd2 (dashed) [upper left panel] or (solid) and 
(dashed) [other three panels] using the polynomial (black), Fourier (gray), and b-spline methods. Top row and lower left panels: 
The 3D spline basis reconstructions correspond to the colored curves, where the purple, blue, orange, and red colors correspond 
to spline basis sets derived from using 10, 14, 22, and 40 ID spline functions in each dimension. Lower right panel: The 2D 
spline basis reconstructions correspond to the colored curves, where the purple, blue, orange, and red colors correspond to 
spline basis sets derived from using 50, 100, 200, and 300 ID spline functions in each dimension. 


10® lO' 10^ 10^ lO'* 



FIG. 3: 2D b-spline fits to (solid) and (dashed), 
where the purple, blue, orange, and red colors correspond to 
spline basis sets derived from using 50, 100, 200, and 300 ID 
spline functions in each dimension. 


can reconstruct the original non-Bunch-Davies shapes. 


The spline expansion method is a new alternative choice 
of basis, that we implement here for the first time, and 
can be used at a reasonable computational cost to obtain 
cosines that can be very close to unity. 

We find that the polynomial basis is good for describ¬ 
ing some non-Bunch-Davies shapes with large features, 
and generally performs better than the Fourier basis. For 
the rest of the shapes we considered, assuming the low 
cosines will steadily increase with the addition of more 
modes, we find that numerical difficulties prevent us from 
generating enough modes to see higher cosines. For most 
shapes, the spline basis performs as well as the polyno¬ 
mial one when equal numbers of modes are chosen. The 
spline basis is both numerically simpler to compute, and 
does not require orthogonalization due to its localized 
nature. This allows a larger number of modes to be ef¬ 
ficiently calculated to improve the match between the 
templates and the actual shapes. 

The spline basis expansion method is very flexible, and 
there are other ways to adapt a spline basis to target spe¬ 
cific shapes better that we have not explored here. For 
example, while our spline bases are derived from equally 






































spaced knots, we note that one can create a basis using 
unequally-spaced knots, such that different fc-space re¬ 
gions are sampled more finely than others. Non-Bunch- 
Davies shapes, where the features are sharply peaked and 
localized near flattened and squeezed configurations, can 
potentially be probed more efficiently by optimizing the 
spline basis in this way. 

Finally, while our analysis takes place in primordial 
/c-space, the flexibility and computational simplicity of 
the spline basis approach translates directly to multipole- 
space, where it complements existing approaches, such 
as the polynomial basis. Utilizing a variety of basis ex¬ 
pansion techniques ensures that the exquisite CMB data 
available now, and in the future, can be used efficiently 
to explore the full theory space of viable inflationary sce¬ 
narios. 
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We can write in terms of its corresponding mode 

functions g^{vi) as 

^fc(^) = ■ (AS) 

For Bunch-Davies modes we choose the positive fre¬ 
quency solution at early times, and using 0 ( 77 ) = 

— l/{r]H) at leading order during inflation, the mode 
functions are given by 


fk iv) 

9ki.g) 


1 


■s^ 


^—iCskrj 



‘■sk^-iCskri 



(A4) 

(AS) 


The choice of Bunch-Davies vacuum can then be ex¬ 
pressed as the following relationship between the field 
and its conjugate momentum in the infinite past. 
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Appendix A: Power spectrum for general initial 
states 

In this appendix we show that the result for the Green’s 
function for general initial states in section |TT] with a 
scale-dependent initial time ?7o(fc) gives the correct form 
of the late-time power spectrum in |10llll| . We start with 
developing a precise definition for the adiabatic vacuum 
at 770 —t —00 or the Bunch-Davies vacuum. Let us write 
the field Xkiv) (defined via ^ terms of an¬ 

nihilation and creation operators at the time 770 , 0 ^( 770 ) 
and ^(^o)) and the mode functions as 


T^kim) = i-iCsk)Xk{vo) ■ (A6) 


Note that this does not imply that the position and mo¬ 
mentum operators commute at all times; it is merely a 
statement of how they are related at 770 —t — 00 . Equiva¬ 
lently, in terms of the mode function ( 77 ) we can write 


d/; 


d77 


V=rio 


1 da 
ad77-''= 


--^fkiVo) = i-ic,k)f^{r]o), (A7) 


for 770 —t — 00 . The above condition is the definition of the 
adiabatic vacuum in the infinite past, or equivalently the 
Bunch-Davies vacuum. As shown in [iniiisi this choice 
corresponds to a minimum uncertainty state. 

The prescription to choose an adiabatic vacuum at a 
finite initial time is to enforce the same condition in eq. 
(A7l at a given 770 . This corresponds to a state which 


minimizes the uncertainty at 77 = 770 . Ghoosing the initial 
density matrix to be unity, but still having a non-Bunch- 
Davies initial state by allowing /?^ 7 ^ 0, i.e. with {rj) 
given by eq. (101 with Ak = 0, this leads to the following 


relation between the Bogoliubov coefficients. 




2cskrio 


^- 2 ic,kr]o 


(AS) 


Gombining this with the usual condition — = 1 

we find that 


Xkiv) = ai:{r]o)f>{r]) + al^{T]o)f<{r]). (Al) 



4c^A:^77g -I-1 


(A9) 


The conjugate momentum, 7i'g(77) can be obtained from 
the quadratic Lagrangian for Xk^v)'! leading order this 
is given by 


Let us now choose 770 to be a function of k such that the 
physical momentum Csk/a{r]Q) crosses some fixed high 
energy scale A of new physics at 770 , then 




dx 


dx% 1 da 

d77 a d77 ’ 


(A2) 


A 

Vo = - 


Hcgk 


(AlO) 
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With the above equations we can write the following so¬ 
lution for and (3^ (note that we are free to choose 
any overall phase), 


a 


> 

k 


P 


> 

k 


\ 2A/i? J 



(All) 

(A12) 


which includes a scale-dependent oscillatory term. The 
argument of the oscillations is usually written as being 
proportional to \n{k/kp), where kp is some fixed pivot 
scale. This can be seen by expanding H around the pivot 
scale so that H[k) « H{kp)[\ — £{N — Np) + ...] « 


H{kp) 


1 -elnf 


, leading to logarithmic oscilla¬ 


tions in the power spectrum. 


For K/H 3> 1 we can now write the late-time power 
spectrum as 


V^k) 


kP Cg 1 
271^ 2e a?{rf) 

— [i- 

STT^ecg 


fkiv)fkiv) 


H 

— Sin 

A 



Cs k—aH 


, (A13) 


Appendix B: Bispectrum for general initial states 


The functions ■k'^^ 2 , and in eqs. 

are given by 


(17l-(19 




1 

2 giCsiCirjo 


2iCsT]o 2 2^ 

kik2k3 

h 

1 

u? 

K 3 


(Bl) 


and 


•^^(;2(pi,P2,P3,?7o) = 


1 


kik2k3 


-2KfKi + KfKl -f KiK^Kl 


gZCaifilJo 


{2KfKl - K^Kl - KiKlKl + A? A|7?o) 


(B2) 


^aacr{Pi,P2,P3,Vo) = ,, , , .3 

[kik 2 k 3 ) L 

giCaifi rjo 


Kf - iKfKl - KfKl + 2KlK^ + 2KiK^Ki 

A? 

+ ^ {k! - 2KlKl) + c, + KlKl - 2KfK^ - 2AiA|A|) r,o 

- tcl {KtKl - 2KlKlKl) 77 ^ 


(B3) 


where for brevity of notation we have suppressed the ex¬ 
plicit momentum dependence of the functions ATi, 7^2, 
and K 3 , 


Kl{pi,P2,P3) 

K2{P1,P2,P3) 

K3{P1,P2,P3) 


Pi+P 2 +P 3 , 

(B4) 

{PlP 2 + P 2 P 3 + P3Piy^^ , 

(B5) 

{PlP2P3y^^ ■ 

(B6) 


1. Flattened limit 


The bj-^ak^CLi;^ term leads to an enhanced flattened 
limit (fci « k 2 + ks, ki being the largest momentum 
mode) bispectrum. Let us first assume that is purely 
imaginary, so that the bispectrum is proportional to the 
imaginary part of 3F. The corresponding T function we 
consider is. 


In the next two subsections we show how one obtains 
the flattened and squeezed enhancements from the func¬ 
tions T (also see jSl]). We work with the simplest func¬ 
tion , though the results are similar for the other two 
functions (or at least for their appropriate sum) as well. 


J'^3(-fci,A:2,fc3,r7o) 


g2Csfci?7o 

+ - 

kikik2k3 



2 

k\kik 2 k 3 
2fCs??o 
fci 



(B7) 
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where fci = —ki + ^2 + k^- In the limit of ki —>■ 0 , the 
exponential can be expanded as 


lim = 1 + iCshVo - — 

fei->-o 2 


6 


(B 8 ) 


Using this in eq. (B7l, and noticing that any term with 
ki in the numerator goes to zero, we find that 

i 


lim J'^ 3 (-fci,fc 2 ,fc 3 , 77 o) = - 

fci->-0 


3(fc2 + k'i)k2k3 


■ 


(B9) 


For fixed 770 , we can set CsIt^qI = 1/A:*. The above limit 
of the three-point function is therefore enhanced (though 
not divergent) in the flattened limit. For r 7 o(fc) with large 
A/H as well we see an enhancement in the flattened limit. 
If we instead assume that is real, then the bispectrum 
is still enhanced, though not in the exactly flattened limit 
but in a near-flattened limit. 


2. Squeezed limit 


Let us now look at the piece. 


*F^3(fci,-fc2,A:3,r7o) = - 

^ ^^Csk2no / 2 2iCsr]o 
k2kik2kz \ k2 


fc|fci/c2A;3 


- cW, , (BIO) 


where k 2 = ki — k 2 + k^. In the squeezed limit {k^ <?; 
ki ~ ^ 2 ) we have k 2 —>■ fcmin) where fcmin is the smallest 
momentum mode observable today. Using this in eq. 
(BIO) we And that 


lim *F^ 3 (fci,-fci, ^ 3 , 770 ) 

^2 ^^min 


2 


piCg fejjjjn77o 
1.2^2 


^min 


2 iCsr]o 

h ■ 

^min 



(Bll) 


For fixed 770 and in the limit of fcmi„ ^ fc*, the 0^770 term 
gives the largest contribution. This term is multiplied 
with a highly oscillatory function though, and averaging 
over the large argument of the cosine (real part of the ex¬ 
ponential) we expect its contribution to vanish. The lead¬ 
ing order contribution is then proportional to 1 / (fc^fc^;„), 
which shows a strong squeezed limit enhancement. In the 
limit of A:min ^ fc* or for 770 (fc) this argument no longer 
holds and we may or may not see enhancements.^*^ 


In our fits with fixed 770 in section [Tv] we took fc* to be similar to 
kmin (~ 10“® Mpc“^) which is what the Planck team had used 
in 1261 . 


Appendix C: B-splines fitting algorithm 


In this appendix, we build on the discussion of sec¬ 
tion |III| to describe in more detail the b-spline fitting 
algorithms we have used, and illustrate b-spline fitting 
examples for a simple ID function and a 2D representa¬ 
tion of the scale-invariant enfolded template. Snippets of 
Mathematica 10 codes we have implemented are shown 
here, and the same algorithms were first constructed and 
explicitly shown as Matlab code in m- 

Our first example considers a spline fit to a simple ID 
function f{x), with 0 < a; < 1. To generate a fit, we first 
make two choices: a choice of basis and a choice of data 
points to fit. The spline basis is determined by a choice 
of fc +1 equidistant knots, {a:o, a;i,..., Xk}, at which each 
basis function’s degree q polynomial pieces will be joined. 
The b-splines in general do not have to be generated using 
constant knot intervals, but for simplicity we always start 
with equidistant knots, and additionally include q extra 
knots at each of a; = 0 and a; = 1 to produce a “clamped” 
basis. Without these extra knots, the generated basis sets 
do not have splines with non-zero amplitudes at a; = 0 
and a; = 1, and it will be difficult to fit functions that are 
non-zero at the endpoints. 

Given these inputs, existing codes, such as Mathemat- 
ica’s BSplineBasis function, can recursively generate a 
basis of k + q b-splines such that each one is spanned by 
q + 2 knots, made up of 5 -|- 1 polynomial pieces, with 
derivatives continuous up to order q — 1. In addition, 
the sum of all b-spline amplitudes at any x is unity. For 
example, the b-splines in fig. are easily generated by 
choosing q = 3 and fc = 3 such that the knots vector 
is knots = {0,0,0,0, 1/3,2/3,1 , 1 , 1 , 1}, and execut¬ 
ing BSplineBasis [{q,knots}, i ,x], where i is an inte¬ 
ger 0<7< k + q — 1 identifying each particular b-spline. 
In this work, we vary the number and widths of the b- 
splines by varying k, but always keep the degree fixed to 
(7 = 3. 

Next, the choice of data points depends on how finely 
we wish to sample f{x). We would like to choose a data 
set consisting of M data pairs, {xi,fi), where ft = f{xi), 
such that our data resolves any potentially fine features 
in the function we would like to fit, without including so 
many extra data points that our numerical calculation 
becomes intractable. The Anal At should ultimately be 
insensitive to the sampling we have chosen. Again, for 
simplicity, we always use equidistant sampling points Xi 
in our analysis, but vary the density of sampling points 
by changing M. 

The spline fit then approximates the original function 
as 


fix) 


N-1 

n—Q 


(Cl) 


where the expansion coefficients {an} are solved for by 
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minimizing the least-squares function, 


M 


LS = Eu.- 


Z=1 


N-1 

E 

n=0 


(^i) 


(C2) 


This requires solving the linear system of equations given 
by B^Bd = B'^ f, where B is wa M x N matrix contain¬ 
ing the values of Bn{xi), and is easily performed with 
algorithms such as Mathematica’s LinearSolve. 

We note one extension of b-splines, called p-splines, 
that aims to avoid overfitting a set of input data by im¬ 
posing smoothness on the resulting fit. Short for “penalty 
b-splines”, in the p-spline method, the fit’s expansion co¬ 
efficients are determined by both the choice of basis and 
the choice of input data, plus a choice of penalty function 
that generally disfavors fits with large differences between 
coefficients of neighboring b-splines . In this context 
one would instead minimize 


M / Af-l \ ^ N-l 

LS = E - E + A ^ , 

i—l \ n—0 / j—k 

(C3) 


where A is a constant that controls the smoothness of the 
fit and k is the order of the penalty, a typical choice being 
fc = 2, such that = aj — 2aj-i + ctj- 2 - The use of 
a penalty is optional, and its main purpose in the con¬ 
text of data fitting is to avoid fitting any noisy features 
in the data. Further, if there are not sufficiently many 
data points sampling f{x), with many more splines than 
data points, then without a penalty the fits may display 
spurious features, as we will deliberately try to show in 
the ID example that follows. 

In fig. 1^ we show fits to f{x) = sin(10a:)/(10a::) 
with different choices of data points and smooth¬ 
ing parameter A. We have fixed the knots at 
{0,0, 0, 0,0.1, 0.2,..., 0.8,0.9,1,1,1,1}, yielding a basis 
of 13 b-splines. We find, as illustrated in the figure, 
that we can achieve good fits without introducing ex¬ 
tra smoothing through a non-zero value of A, as long as 
we fit to enough data points. So we now continue to an 
example of fitting a primordial shape in two dimensions, 
without a penalty. 

To illustrate the b-spline fitting algorithm in two di¬ 
mensions, we construct a basis of 2D b-splines and use it 
to fit the scale-invariant enfolded template, 

Senf{x,y) = -{l-x-y-x'^-y'^+x^+y^ 
xy 

- x^y - xy"^ + 3xy ), (C4) 

where x = k^/ki and y = k 2 /ki. As in the ID case, the 
inputs to the fitting algorithms are made up of a choice 
of basis and a set of data points. The basis is specified by 
a choice of polynomial degree and a sequence of knots in 
each of the two dimensions, x and y. In our particular ap¬ 
plication, since we are aiming to fit shape functions that 



0.0 0.2 0.4 0.6 0.8 1.0 


X 


FIG. 4: Example ID b-spline fits to f{x) = sin(10a;)/(10a;). 
All fits shown have used a basis of 13 splines. In the up¬ 
per panel, the fits have been computed using 6 data points, 
sampled at Xi = {0, 0.2, 0.4, 0.6,0.8,1}. Attempting to con¬ 
struct such a fit without smoothing [solid red] causes the 
splines to produce spurious features, especially at smaller val¬ 
ues of X, while introducing a penalty and a small amount 
of smoothing [dashed red], A = 0.1, restores the fit to a 
reasonable representation of the true function. The lower 
panel is the same fit, except with 21 data points, sampled 
at Xi = {0,0.05, 0.1,..., 0.9,0.95,1}. In this case, the fits 
with [dashed blue] and without [solid blue] smoothing are 
very similar. 


are symmetric in their wavenumber arguments, we only 
specify the knots and degree in one dimension, and use 
the same b-spline basis for the additional second dimen¬ 
sion. The 2D b-splines are then made of products of any 
two ID splines, for example, Bn{x)Bm{y)- The data are 
given by {xi,yj, S{xi,yj)) and stored in Ytj = S{xi,yj), 
where again some care must be taken in the choice of 
sampling, which must be dense enough to capture any 
small features such as oscillations that we would like to 
capture in the resulting fit. 

The 2D analogue of the least-squares function in eq. 


(C2l is 


MM/ N-l N-l 

LS = ^ ^ ^ ^ I Yij ^ ^ ^ ^ 0^mnBm{Xi)Bn{yj 

i—l j — 1 \ m—0 n—0 


(C5) 


To turn the problem of solving for the expansion coeffi¬ 
cients Omn into a linear system, we create a regression 
basis C from the M x N b-spline basis matrices in each 
dimension, Bi and i? 2 , which in our case are equal. We 
define 


C = (Ba 0 ej) © (e^ 0 Bi) = SiDBa , (C6) 
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where 0 is the Kronecker product, © is an element-by¬ 
element multiplication, the second eqnality defines the □ 
operation, and and ck are vectors of I’s with length 
L each. In Mathematica, we define the □ operation as 
box: 

box[Bl_,B2_] :=Module[{K,L,eK,eL} , 

K = Length[B1 [[1 , All]] ] ; 

L = Length [B2 [ [1,All]]] ; 
eK=ConstantArray[1,{K}]; 
eL=ConstantArray[1,{L}]; 
KroneckerProduct[B2,{eL}] 
*KroneckerProduct [{eK},Bl] 

] 


0 10 20 30 40 50 


oq. 




o 

o 


0 10 20 30 40 50 

Number of modes 



Then by stacking the columns of the coefficients array 
amn and the data array Yij to get vectors /3 and y respec¬ 
tively, the task of finding a solution for the coefficients 
is once again rednced to solving a linear system of eqna- 
tions given by IT(7/3 = C"'"Wy. Here IT is a matrix 
containing weights, which may be different for each data 
point, but for simplicity we restrict ourselves to using a 
weight of unity for all of our data. 

For modest amounts of data and numbers of b-splines, 
one can qnickly solve for the coefficients in this straight¬ 
forward way. For large data sets and nnmbers of b- 
splines, however, this approach becomes computationally 
cnmbersome due to the large size of C. While it is still 
possible to numerically solve for the coefficients amn us¬ 
ing a low-level language like C(++) or Fortran, we have 
instead nsed algorithms developed for higher level lan¬ 
guages, such as Matlab in (51], nsing only vector and 
matrix operations. This has the benefit of being easier 
to implement, while still being able to sidestep much of 
the memory storage and speed issnes typical of a more 
brnte-force approach in a high-level langnage. Instead 
of starting with a calcnlation of C in the brute-force ap¬ 
proach, the algorithm from m that we have adopted 
computes (7^ IT(7 and C"^Wy using only B. We refer the 
reader to m for a detailed discnssion of how the method 
itself is devised and constrncted, or to see the equivalent 
Matlab code, and present here an implementation of the 
b-spline fitting algorithms in Mathematica. 

The normal equations can be efficiently constructed 
and solved, given an inpnt of data in Y and information 
abont the data sampling and b-spline basis in B: 

get2dfit [Y_,B_] :=Module [{m,n,W,R,r, 

F,a,A}, 

m=Length[Y[[l,A11]]]; 
n=Length[B[[l,A11]]]; 

W=ConstantArray [1; 

R=Transpose[B].(W*Y).B; 
r=ArrayReshape[R,{n*n,l}]; 

F=Transpose[box[B,B]].W.box[B.B]; 
F=ArrayReshape[F,{n,n,n,n}]; 

F = TensorTranspose[F,Cycles{{3,2}}]] ; 

F = Arr ay Reshape [F,-[n*n,n*n}] ; 
a=LinearSolve[F,r]; 

A=ArrayReshape[a,{n,n}] 


FIG. 5: Cosines between the enfolded shape and the 2D b- 
spline fits generated by a basis with 10 splines per dimension. 

] 

After the matrix of coefficients is solved for, we con¬ 
struct the final fit through two steps. First, we map 
the coefficients output as A from get2dfit to a new 
set of coefficients that corresponds to a 2D basis of 
splines which is symmetric in its two arguments, so 
that each 2D basis mode is a snm of np to two terms: 
Bi{x)Bj{y) + Bj{x)Bi{y). Second, in building np the fit, 
mode by mode, we start with the modes that contribnte 
most to the fit. Since the b-splines in a choice of basis 
have similar shapes and amplitudes, we use the magni¬ 
tude of the Aij coefficient as a proxy for gauging how 
much any particnlar b-spline contributes to a fit’s overall 
cosine with the original shape. This motivates building 
np a fit by adding in modes, starting with those that 
have the largest \Aij\. The cosine is then computed in 
the nsnal way, throngh an inner prodnct over (x, ?/)-space 
between the original shape S'enf and the fit. 

As an example, we nse a basis of 10 splines in each 
dimension constructed by choosing k = 7 and g = 3, 
and nse as onr data set a grid of uniformly spaced (x, y) 
valnes from taking 50 samples in each dimension, to com¬ 
pute the matrix of coefficients Aij, nsing the algorithms 
box and get2dfit above. The total number of symmet¬ 
ric modes is then 55, and the modes are ordered by their 
corresponding values of largest to smallest | Aij \ to pro¬ 
duce the cosines in fig. A visual comparison of the full 
fit using 55 modes and the original enfolded template is 
given in fig. 

For our 3D fits, the fitting method is the same: a choice 
of b-spline basis and data set make up the inputs to the 
algorithm, which returns an array Amnp containing the 
expansion coefficients that approximate the input shape 
as S {k\, k 2 , ki^') — Ammn Bnr}.{k-\ ^Bn,(k^^BnjkA. 

However, due to the higher dimensionality of the prob¬ 
lem, we mnst introduce a new fnnction, rho, to gen¬ 
eralize the matrix product to the product of a matrix 
and a 3D array. Below, we show Mathematica code for 
rho [A, B, p], which computes the normal matrix product 
between rows of A and the p**' column of B, resnlting in 
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FIG. 6: The original template Seni{x,y) [upper panel] and 
its b-spline fit using a set of 55 modes [lower panel]. The 
amplitudes range from 0 (violet) to 1 (red) in both panels. 

a product, C, which has the same dimensions as B: 

rho[A_,B_,p_] :=Module [{sa,sb,n,ip, 
cycles ,sbip,prodsbip ,tempB ,C}, 
sa=Dimensions[A]; 
sb=Dimensions[B]; 
n = Length [sb] ; 

ip=Join[Range[p+1,n],Range[1,p-1]]; 
Which[ 

p = = l,cycles = Cycles [{}] , 
p==2,cycles=Cycles[{{l,3,2}}], 


p==3,cycles = Cycles [{{1 ,2,3}}]] ; 
tempB=TensorTranspose[B,cycles]; 
sbip = sb [ [ip]]; 
prodsbip = Product[ship [[i]] , 

{i,1,Length[ship]}] ; 
tempB = ArrayReshape[tempB,{sb[[p]] , 
prodsbip}]; 

C=Transpose[A].tempB; 

C=ArrayReshape[C,Join[{sa[[2]]}, 
sb [[ip]]]] ; 

C = TensorTranspose [C , 

InversePermutation [cycles]] 

] 

Given this definition of rho, the 3D b-spline fit coeffi¬ 
cients are calculated using get3dfit: 

get3dfit [Y_,B_] ;=Module [{m,n,W,F,R,A}, 
m=Length[Y[[1,1,All]]]; 
n=Length[B[[l,A11]]]; 

W=ConstantArray[l,{m,m,m}]; 

F = rho[box [B,B] ,W , 1] ; 

F = rho[box[B ,B] ,F,2] ; 

F = rho[box[B,B] ,F ,3] ; 

R=rho[B,Y*W , 1] ; 

R=rho[B,R , 2] ; 

R=rho[B,R , 3] ; 

F = Arr ay Reshape [F,-[n,n,n,n,n,n}] ; 
F=TensorTranspose[F, 

Cycles [{{3,2,4,5}}]] ; 

F=ArrayReshape[F,{n"3,n“3}]; 

A=LinearSolve[F, 

ArrayReshape [R,{n“3,1}]] ; 

A = Arr ay Reshape [A,-[n,n,n}] 

] 
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